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Abstract. Usual approach to the foundations of quantum statistical physics is based 
on conventional micro-canonical ensemble as a starting point for deriving Boltzmann- 
Gibbs (BG) equilibrium. It leaves, however, a number of conceptual and practical 
questions unanswered. Here we discuss these questions, thereby motivating the study 
of a natural alternative known as Quantum Micro-Canonical (QMC) ensemble. We 
present a detailed numerical study of the properties of the QMC ensemble for finite 
quantum systems revealing a good agreement with the existing analytical results for 
large quantum systems. We also propose the way to introduce analytical corrections 
accounting for finite-size effects. With the above corrections, the agreement between 
the analytical and the numerical results becomes very accurate. The QMC ensemble 
leads to an unconventional kind of equilibrium, which may be realizable after strong 
perturbations in small isolated quantum systems having large number of levels. We 
demonstrate that the variance of energy fluctuations can be used to discriminate the 
QMC equilibrium from the BG equilibrium. We further suggest that the reason, why 
BG equilibrium commonly occurs in nature rather than the QMC-type equilibrium, 
has something to do with the notion of quantum collapse. 



1. Introduction 

This article is about the properties of a quantum statistical ensemble, which is, in a sense, 
opposite to the conventional micro-canonical ensemble. It is an extension of the previous 
work of one of uspQ. The present work contains a discussion of underlying conceptual 
issues and detailed analytical and numerical investigations of finite quantum systems. 
These two aspects of the article are complementary but otherwise quite different in 
scope and style. Sections [2] and [8] are dealing with conceptual issues. They motivate 
and outline a broader agenda. The remaining sections flHKplus Appendices) are more 
specialized. We hope that the conceptual part would motivate the readers unfamiliar 
with the above agenda to look into the rest of the article. 
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2. Narrow eigenstate energy window in conventional statistical ensembles 
and the motivation to look beyond it 

The standard quantum-mechanical derivation of Boltzmann-Gibbs (BG) equilibrium is 
based on the conventional micro-canonical postulate formulated as follows. One should 
consider a subsystem in a macroscopic environment, assume, that the two form an 
isolated macroscopic system, neglect the interaction between them, and then select 
the statistical ensemble amounting to an incoherent mixture of eigenstates of the total 
system in a narrow energy window around a given value of energy, E av . The width 
of this micro-canonical energy window A mc should be large enough to include a very 
large number of quantum levels, but, at the same time, small enough to imply the 
negligible change of temperature across this energy window. This condition can be 
formally expressed as follows: 

A - >> Rt)' (1) 

and 
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Here v(E) is the density of energy states of the isolated macroscopic system, which 
should be approximately equal to the density of states of the environment at the same 
energy, z/ is unimportant normalization constant, and T is temperature in energy units. 
Given the exponentially large number of levels of a macroscopic system, the above two 
inequalities can be easily satisfied for A mc ranging over many orders of magnitude. 

Recently, further progress was made by showing that a single pure state of a 
macroscopic system formed by a coherent superposition of eigenstates randomly selected 
from the above energy window would also lead to a density matrix of the subsystem 
exhibiting the BG distribution [2j El HJ [5]. It has been further conjectured that, if the 
isolated macroscopic system is chaotic in the classical limit, then even a single eigenstate 
would lead for the BG distribution in the density matrix of the subsystem[6]. 

Yet the question remains: What is the justification of the narrow energy window 
in the microcanonical postulate? The origin, the limits and the implications of this 
condition are not well understood at present. 

The microcanonical postulate as such remains a rather formal construction, because 
in reality we do not encounter macroscopic systems isolated from the rest of the world. 
Even, if we were to find one, the non-relativistic quantum-mechanical description of such 
a system normally leads to unrealistically small spacings between adjacent eigenstates 
- much smaller than the minimal natural width of an energy level estimated as fr/tjj, 
where ty « 4 X 10 17 s is the Big-Bang age of the Universe. The real width of the energy 
levels should, of course, be much larger due to relativistic effects. Although the notion 
of an energy eigenstate of a macroscopic system is not well-defined, one can still hope 
that the statistical treatment based on this notion would produce an outcome, whose 
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validity extends beyond the systems with well-defined energy eigenstates — in the same 
way as BG statistics derivable for classical systems extends to quantum systems. 

Considering systems with well-defined eigenstates, one can proceed by giving the 
conventional microcanonical ensemble with its narrow energy window the status of a 
fundamental postulate. However, such an approach to justifying BG equilibrium is not 
quite satisfactory. If instead of simply postulating the BG equilibrium for quantum 
systems, one were to claim its derivation from a more fundamental postulate, then the 
introduction of this extra postulate would be justified only if, in addition to the BG 
equilibrium, it would lead to an explanation of some other previously unexplained facts 
about nature and predict new ones. The micro-canonical postulate, however, does not 
lead to any other observable prediction besides BG equilibrium itself. 

It could have been that the narrow energy condition is only a simplification, while 
the BG distribution emerges for a broader class of ensembles. However, as discussed 
in Refs.[7J [TJ [8], the departure from the narrow energy condition easily spoils the BG 
statistics for a subsystem. It leads to a mixture of different thermal distributions, which 
does not produce a distribution characterized by a single temperature. Yet one could 
have asked the question: If it is not a narrow energy window, then what is the most 
likely shape of the participation function for the energy eigenstates of the entire system? 
Can it be that, if the most likely shape is found, then the resulting statistics for the 
subsystems is of BG type? 

One can, in particular, recall that the conventional canonical ensemble does not 
appear to have narrow energy window for participating eigenstates. To each eigenstate 
with energy E, it assigns probability weight p(E) = exp(—E/T ), which follows either 
from maximizing entropy — ^^(-E^log p(Ei) subject to the total energy constraint, 
or from the practical consideration that, if one slowly isolates a large quantum system 
from its environment, which initially was in BG equilibrium with temperature T , then 
the diagonal elements of the resulting density matrix will have the above form. Here, it 
should be noted, however, that the quantity of practical interest in not the probability 
weight function p(E) as such but rather its product with the density of states p(E)v(E). 
Since p(E) decays exponentially with constant "rate" 1/T , while for a macroscopic 
system v{E) increases exponentially with decreasing "rate" 1/T(E), where T(E) is 
defined by Eq.(|2]), the product would have a narrow maximum at energy E = E- 
corresponding to T(E av ) = T . The width of that maximum A c is determined by 
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where N p is the number of particles in the system. This gives A c ~ To^/Wp. As a result, 
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i.e. condition ([2]) is satisfied. Thus the canonical ensemble for an isolated macroscopic 
quantum system, in effect, imposes a narrow energy window condition. 
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The idea of maximization entropy at a given energy is closely related to the notion 
of chaos and ergodicity in the phase space of classical systems. For quantum systems, the 
notion of phase space trajectories is not defined because of the Heisenberg uncertainty 
principle. Hence the maximization of conventional entropy cannot be justified on the 
basis of the dynamics of isolated quantum systems. 

To address the issue of what may happen beyond the narrow energy window 
imposed by the conventional ensembles, an ensemble with unrestricted participation of 
eigenstates and without the requirement of the maximization of conventional entropy, 
referred to below as Quantum Micro-Canonical ensemble [9j [TOj El EE], is a conceptually 
necessary limit to consider. (Similar quantum canonical ensembles have also been 
defined Refs.0 [Lj.) 

In addition to the above conceptual issues, the basic concern in the context of 
the foundations of quantum statistical physics remains about the actual behavior of 
an isolated quantum system perturbed from equilibrium and not interacting with a 
large thermal bath. It is overwhelming intuition of many researchers, including the 
present authors, that an isolated macroscopic classical system with generic non-linear 
interaction between particles would evolve dynamically to exhibit the BG statistics 
for small subsystems. The situation is, however, not so trivial for quantum systems, 
because they can be in a superposition of states having different total energies. If, after 
the perturbation, an isolated quantum system is in a broad superposition of eigenstates 
violating the narrow energy condition, then the occupations numbers of eigenstates will 
not change with time under the dynamics guided by linear Schrodinger equation, and 
hence the narrow energy condition will never emerge. 



3. Quantum micro-canonical ensemble 

Generic wave function \I/ of a quantum system is a superposition of eigenstates (ftf. 

N 
i=l 

where are complex amplitudes, and N the total number of eigenstates. For an isolated 
quantum system, fixing energy E av associated with wave function \1/ implies 

N 

£ E iVi = E av , (6) 

i=l 

where pi = \Ci\ 2 are the occupation numbers of quantum eigenstates, and are the 
corresponding energies. We shall refer to E av as "average energy" [not to be confused 
with the uniform average of all Ej\. 

In this work, we focus on Quantum Micro- Canonical (QMC) ensemble [9] fTU | ITT | [T] . 
where all quantum wave functions satisfying average energy constraint fl6]) are equally 
probable in the Hilbert space of the system. The QMC ensemble amounts to a natural 
generalization of the equipartition postulate for energy shells in classical phase spaces. 
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Despite this analogy, the QMC ensemble is quite different from the conventional micro- 
canonical ensemble, because it does not limit the participating eigenstates to a narrow 
energy window around E &v . 

It was shown in Ref. [1] that the QMC ensemble does not lead to the BG statistics, 
at least in the absence of interaction between the subsystem and the environment. This 
negative result is important, because it forces one to think seriously about justifying 
the narrow energy window in the conventional microcanonical ensemble. As discussed 
in Refs. [U |8] and in Section [BJ quantum collapse may play here a fundamental role. 
Further study of the QMC ensemble is also important, because this ensemble may be 
realizable in small isolated quantum systems with a large number of quantum levels. 

Beyond the subject of foundations of quantum statistical physics, the properties 
of the QMC ensemble are also relevant to the studies of various properties of Hilbert 
spaces, e.g. in the context of quantum computing or many-particle entanglement. 

It is finally plausible on general grounds that such a natural quantum ensemble 
would find many realizations not discovered so far. For example, it was demonstrated 
in Ref.p] that the expansion of the basis states of sparce random matrices in terms of 
eigenstates exhibits the QMC type of statistics. 

The description of the QMC ensemble involves two steps: (1) characterization of 
the shape of the energy window function — in this case, the participation function for 
the eigenstates of the entire isolated system; and (2) obtaining the density matrix for 
the subsystem of interest. Sections HH6] below deal with step (1), while Section [7] deals 
with step (2). 

In Section HI we summarize the analytical results obtained in Ref. [1] for the QMC 
ensemble in large- iV isolated quantum systems. In Section [5J we propose a method of 
introducing finite- iV corrections to the large- N approximation. In Section [61 we present 
the numerical tests of the resulting description for a 12-level isolated quantum system. 
In Section [7J we examine numerically the density matrix for a subsystem of a 12-level 
system, and compare the QMC-based result with that of the conventional canonical 



ensemble. Section [HJ discusses open questions motivated by our results. |Appendix A 



describes a practical numerical algorithm for finding an important characteristic function 
appearing in the course of the analytical approximation procedure. Appendix B gives 
the details of our Monte-Carlo sampling algorithm for the QMC ensemble. 



4. Description of the previous QMC results 

We consider quantum spectum {Ei} consisting of iV energy levels. The values of energies 
are assumed to be ordered with Emm = E\ and Emax = En- 

The analytical approach and the analytical results of Ref. pQ can be described as 
follows. 

The QMC ensemble is defined by specifying uniform joint probability distribution 
for complex amplitudes = \Ci\e %ipi in the Hilbert space — subject to average energy 
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constraint OH]) and to the normalization constraint 

JV 

5> = 1 " (7) 

i 

The two constraints are of the second order in variables |Cj| and do not depend on (pi. 

To avoid dealing with curved manifolds, we now change from variables |Cj| to 
variables Pi = \Ci\ 2 . Completely random choice of states in Hilbert space implies that 
the probability measure is proportional to the volume element: 

N 1 N N 

dv^HiQi diai = ~Y[dtpi ddai 2 ) = ~]Jdw d Pt . (8) 

i i i 

Since constraints (EJ U} do not depend on ip i} the assignment of variables ipi is totally 
random. The resulting uniform distribution of (pi can be just factorized out and then 
integrated. What remains is the Euclidean subspace of variables pi with uniform joint 
probability measure on the manifold constrained by conditions ([6], [7]) together with the 
positivity condition 

Pl > 0, VI (9) 

In the pj-space, due to the linear nature of constraints (l6|7|9|) . the boundaries of the 
resulting manifold are flat, i.e. the manifold itself is an (iV — 2)-dimensional polyhedron. 
We will refer to it as "QMC polyhedron". Our goal is to obtain marginal probability 
distribution Pk{pk) for a given occupation number p^ and the corresponding average 
value (pk) representing the participation of the kth eigenstate in the QMC ensemble. 

Probability distribution Pk{pk) is proportional to the (N — 3) dimensional volume 
of the intersection of the QMC polyhedron with hyperplane pk = const. (Here we use 
variable pk as a label of the coordinate axis, while in Pk{pk) and everywhere below, pk 
refers to the value of const.) The intersection manifold is to be denoted as Aik, and its 
volume as Vk(pk)- This volume plays the role of unnormalized probability distribution 
such that 

P k (p k ) = i VkiPk) , (10) 

and the average value of pk is 

In the space of all variables excluding pk — to be denoted as {pi}k, manifold Aik 
is defined by conditions: 

N 

jCPi = i-p*, ( 12 ) 

and 

N 

(Ei - E aw ) Pi = -{E k - E av )p k , (13) 
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together with positivity conditions (J5J). 

The most relevant limit to consider now is p k <C 1, because, for N ^> 1, it is 
improbable that a single occupation number p k in a randomly chosen superposition 
of all quantum states becomes comparable to 1. We refer to this limit as "small-p^ 
approximation." It gives pQ 

V k ( Pk ) = V k (0)e' N ^ 1+x ^- E ^. (14) 

where X k is the volume renormalization parameter associated with the shift of energy 
hyperplane ( fi~3i) . 

It was further shown in Ref. [1] that, in the limit iV ^> 1 all parameters A^ are 
approximately equal to the same value — to be denoted as A, i.e. in this limit 

V k ( Pk ) = V k (0)e- Np ^ 1+x ^- E ^. (15) 

The corresponding average value of the occupation number p k is 

« = iv[i+A ( k -jyr (16) 

The value of parameter A can now be found numerically by substituting Eq. flTB)) into 
either averaged normalization condition 

N 

£<P*> = 1, (17) 

k=l 

or averaged energy condition 

N 

J2(E k -E av )(p k ) = 0. (18) 

k=l 

The two values of A are guaranteed to be the same for any N, probably because 
formula ffTB]) is analytically accurate leading term in the large- N expansion. [For finite- 
iV systems, the actual behavior of V k {p k ) can exhibit noticeable deviations from the 
exponential dependence of Eq. (TT3|) ]. 

Behind the geometry-based discussion one should not overlook that Eq. ffl5l) 
represents grand-canonical Gibbs ensemble defined in the Hilbert space, or the space of 
occupation numbers p k . The problem is, in fact, equivalent to finding the probability 
distributions for a system of iV classical states characterized by one-particle energies 
E k , with p k being the number of classical particles occupying the k-th state. In such a 
problem both the total energy and the total number of classical particles is conserved, 
which leads to the grand canonical ensemble in the Hilbert space. This analogy, however, 
does not imply conventional observable results, because it pertains to the probability 
distribution P k (p k ), while, in the original quantum problem, p k itself is the probability 
of observing the system in the kth state, i.e. the prediction of the observable outcome 
requires finding the average value (p k ). It is at this last step that the departure from 
the conventional statistics occurs. In a broader context, this kind of result associated 
with double- averaging is known as superstatistics[13j. 
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The grand-canonical analogy becomes evident, once formula (IT5"j) is rewritten as 
V k (p k ) = V k (0)exp{-NX(E k p k - E xPk )}, (19) 

where 

E x = £av - y (20) 

Equation ([19]) implies that XN is the inverse temperature in the Hilbert space, and E\ 
is the chemical potential. 

The value of A is zero, when E av is equal to E av0 , the uniform average of all energies 
in spectrum {K}[TT1 IT]: 

1 N 

-EavO = ^^£i- (21) 
i=l 

We also recall that A > for E av < E av0 , and A < for S av > E av0 . Here we follow 
the convention of Ref. [TJ, that, unless specified otherwise, the origin of the energy axis 
is set at E av0 , i.e., E av0 = 0. 

The chemical potential E\ is simultaneously the pole of function (p k )(E av ) given by 
Eq. (jT6l) . which implies that for the small-p^ approximation to be valid for all levels, E\ 
should be sufficiently below the lowest energy level for positive A or above the highest 
energy level for negative A. In the former case, "sufficiently below" means thatJJJ 

E min - E\ » ^ av ~ Emin . (22) 

It is shown in Ref.pQ that the above condition is violated for a typical macroscopic 
system at a typical value of E av , which leads to condensation into the lowest energy 
state, i.e. (pi) becomes comparable to one. Condition fT22|) can also be easily violated 
for finite-N systems leading to the departure from the small-p^ results (I15|ll6p for several 
lowest levels. In the both cases, it is necessary to consider what happens beyond the 
small-pfc approximation. 

When typical values of p k are not small, it is necessary to recall [TJ that parameters 
Afc are defined as functions of average energy E av as follows: 

Afc[i?av] = TT , (23) 

where [in this equation only] V k is the volume of the manifold constrained in space {pi} k 
by conditions 

N 
N 

J2(Ei-E av ) Pl = v, (25) 

together with positivity conditions (Q. Here v is the energy shift parameter. We follow 
the convention of Ref.pQ to use square brackets in functions \ k [E av ] and later in A[i? av ]. 
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Given the above definition of Ajj.[£J av ], the volume of manifold M.k defined by 
Eqs. fl9lfT2][T3]) can be expressed as 



V k {p k ) = V*(0)exp { {N - 3) 



p (-Efc--Eav)pfc 
Eav —k 



log(l-p fe )+ / Xk[E]dE 



.(26) 



[This equation does not appear explicitly in Ref.pQ, but otherwise it is an obvious 
intermediate step in the calculation described there.] In the limit N ^> 1, all Xk[E] are 
approximately equal to a single function X[E] defined in interval [E m i U , E max \. Therefore, 
Eq. fT26l) becomes [1] 



V k (p k ) = V k (0)exp{(N-3) 



p (-Efc--Eav)Pfc 
Eav ^ 



log(l-p fe )+ / X[E]dE 



.(27) 



Functions Vk(pk) given by both Eq. (l26|) and (!27j) are not necessarily defined in the entire 
range < p k < 1, because, for too large pk, and, in the case of k = 1 or k = N, for 
too small pk, the intersection manifold defined by Eqs. fl9lfT2"l[T3"]) may simply not exist. 
The practical rule is that Vkipk) is defined, when \k[E] or X[E] are defined at the upper 
integration limit i? av — ( Ek ~Ea.v)Pk . This limit should fall within interval [E m - m , £ max ] 
for Eq. (l27|) and within sometimes different boundaries (given in Section [5]) for Eq. (l26|) . 
The upper cutoffs for pk are given in Appendix D of Ref.pQ and the lower cutoffs in 
Section [5j The mean-A approximation represented by Eq.(T57|) neglects the existence of 
lower cutoffs of pk for k = 1 and k = N. These cutoffs, however, can be incorporated, at 
the level of finite- N corrections based on Eq.( l26l) . If a lower cutoff is present, then V&(0) 
is not defined, and the prefactor in Eq.( !26|) should be viewed as simply a constant, which 
is then canceled by normalization in Eq.f llOp . Also, when the lower cutoff is present, the 
lower integration limit ll E &v n in Eq.f l26p should be replaced with the minimum or the 
maximum value of energy, where Xk[E] is defined, for k = 1 and k = N, respectively. 
In general, function Vk(pk) becomes exponentially small as it approaches either of its 
cutoffs. Therefore, as N increases, the practical significance of these cutoffs decreases. 

In Eq.(|271). we took, at first sight inconsistently, the large- N limit Xk[E] w X[E] but 
did not approximate iV — 3 ~ N. The reason is that we are aiming at making numerical 
tests for Hilbert space with iV ~ 10, where condition iV — 3 ^> 1 is not well fulfilled, 
and, since the term (N — 3)log(l — Pk) is an exact onep], we thereby limit the finite- iV 
error only to the integral in Eq. (|27l) and deal with this error later in Section [5j It is also 
worth noting that the small-p^ formulas f fl4|) and f JT5|) can be obtained by linearizing 
exponents in Eqs. (I26j) and fT2T|) with respect to pk, and approximating iV — 3 ~ N. 

Equation ( F271) was the basis for proving the condensation of the QMC ensemble 
into the lowest energy state of macroscopic systemspQ. It indicated that the character 
of function Vi(pi) was not exponential but rather narrowly peaked around p\ = (pi). 
This condensation is different from the Bose-Einstein condensation, in the sense that 
it occurs in the many-particle ground state independently of the type of the particles 
composing the system. 
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5. Ansatz for deviations of \k[E av ] from A[J5 av ] 

While we cannot propose a controllable procedure for finding corrections for Afc[i? av ] 
with respect to A[i? av ], we introduce here an ansatz, which makes these corrections 
in a conceptually meaningful way and, at the same time, significantly improves the 
agreement with Monte-Carlo results. 

The ansatz is based on the following considerations. Function A[i? av ] is defined in 
interval [E m i D , E max \. It has two universal properties in the large- N limit pQ: 

(i) A[£ av0 ] = 0; 

(ii) for E av < E av0 , A[-Eav] quickly approaches the asymptotic form: 

X[E av ] « 1 , (28) 



y av -"-'mm 



while for E av > E av0 , the asymptotic form is 

X[E av ] » 1 ■ (29) 

-'-'av ^max 

In other words, the behavior of the entire function Af-E^v] is almost entirely 
determined just by the values of three parameters E av0 , E min and E max . 

Now, function A[i? av ] characterizes the entire quantum spectrum {E^}, while 
functions X k [E av ] characterize spectra {Ei} k obtained from {Ei} by removing energy 
level E k . As a result, the uniform average of all energies for spectrum {E{\ k is shifted 
to 

E av ok = E av0 H — — - — , (30) 

while the end points of the spectrum remain unchanged, unless k = 1 or k = N, when, 
respectively, the minimum or the maximum energy levels are removed: 

Emink={ f ; k 1 (31) 

. E max : k < N, 
hnuvsL " Efif-i : k = N. {6Z) 

Apart from the above specific differences Xk[E av ] should be reminiscent of A[£/ av ]. 
Therefore, in order to approximate \ k [E av ] from the knowledge of A[£? av ] we adopt the 
following two-step procedure: 

(I) Shift of X[E av ] to move E av0 to E av0k : 

AJ-^av] = X[E av + E av0 — E av0k \. (33) 

This also shifts E min and E max (the end points of X[E av ]) to E' mink = E min + E av0 - E av0k 
and E' maxk = E max + E av0 - E av0k , respectively 

(II) Uniform contraction or expansion of energy differences above and below 
E av0k to shift E' mink and £^ axfc to, respectively E mink and E maxk given by Eqs.fl3TT|32D 
and simultaneously rescale \' k itself to make sure that, if A' fc [i? av ] has asymptotic 
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forms -b — kr, — above E'- h and ^ — L below E',, then the resulting function 

Esv ~ minfc S av--B maxfe max *' 

Afc[-E av ] has asymptotic forms gav _ff — above E m[n k and g av _]; maxfc below E majX k- This 
transformation has form 



\[E a 

where 



^^fct-^avOfc + &-{E av — -EavOfc)] 
^-A'^fi^avOfc + Ot + (E av — E av Qk)} 



E av < -EavOfc) 
E av > E av Qk, 



(34) 



«- = — p > ( 35 ) 

-^minfc -^avOfc 

a+ = ~p> Z~/r • ( 36 ) 

^maxfc -^avOfc 



and 



The resulting function A&[£a V ] is then used in Eq. fl26|) . 

The last step automatically corrects the drawback of the mean-A formula ( 1271) 
mentioned in Section HI Namely, the mean-A approximation does not have lower cutoffs 
for Vi(p\) and Vn(pn), while the exact solution should have them. Indeed, if _E av falls 
between E m ± n and E 2 , then Vi(pi) has the lower cuttoff 

E2 — E av 

Pimm = p = , (37) 

which corresponds to only two lowest levels E min and E 2 being occupied with the 
occupation numbers determined by the value of E av . Once higher energy levels become 
occupied, the value of p\ becomes greater than pi mm to fulfill average energy condition 
(E]). Likewise, if E av falls between E N _i and E max , then lower cutoff for Vn(pn) is 

En-i - E av 

PNmin = -= = • (38) 

In our approximation procedure, Vkipk) is defined only, when Xk[E] is defined at the 
upper integration limit in Eq. (l26|) . i.e. when E av — ( Ek ~E™)pk falls in the interval between 
E mink and E maxk . The appearance of E minl ^ E min and E maxN ^ E max in Eqs.((35Jl36J) 
then implies the correct lower cuttoff values pi mm and pNmm, respectively. These lower 
cutoffs are unimportant in the large- N limit, but become important for finite- N systems, 
like the one considered numerically in this work. 

6. Numerical investigation of the QMC ensemble for a finite isolated 
quantum system 

Here we compare analytical and numerical results for a quantum spectrum that has 
N = 12 energy levels. This is nearly the maximum number of quantum levels, whose 
statistics we can access in reasonable time with the Monte-Carlo algorithm described 
in Appendix B The values of the spectrum energies are { — 1.27, —0.93, —0.68, —0.47, 
-0.27, -0.09, 0.09, 0.27, 0.47, 0.68, 0.93, 1.27} (also shown in Fig. 1). This spectrum 
is obtained as a strongly discretized version of the Gaussian density of states with 




Figure 1. Main figure: Function A[-E7 av ] in various approximations for N = 12 
quantum spectrum represented by vertical lines. Long-dash magenta line — small-p^ 
approximation based on Eqs. (|15ll6[) ; short-dash green line — mean- A approximation 
based on Eq. (|2"T)k thick solid red line — finite- AT corrected calculation based on Eq. (|2l)|) 
and the ansatz of Section O thin solid black line — asymptotic behavior in the large- 
N limit given by Eqs. (|28l29[) . Inset: Normalization test of the mean- A approximation 
(dashed green line) and finite- N corrected calculation (solid red line). 



maximum at E av o = 0. (Nearly Gaussian density of states is genetically expected for 
macroscopic systems pp.) 

Figure 1 presents curves for A[J5 av ] obtained in various theoretical approximations: 
namely, small-p^ approximation based on Eqs. (jl5)16j) ; mean-A approximation based on 
Eq. (I27j) ; and the calculation with finite- iV corrections based on Eq. (126|) and the ansatz of 
Section [5j The numerical algorithms for computing the latter two curves are described in 



Appendix A Figure 1 also includes curves representing the large- N asymptotic behavior 



(12811291). 

The calculations of all three curves are based on satisfying energy constraint 
(118p . The small-pfc approximation, guarantees that the normalization constraint (jl7jl is 
satisfied at the same time. However, the other two approximations do not guarantee it a 
priori. Therefore, the check of how the normalization constraint is satisfied constitutes 
an indication of the accuracy of the approximations. This check is presented in the 
inset of Fig. 1 as a function of E av . It indicates that the normalization for the mean-A 
approximation falls short of 1 by about 10 per cent (roughly 1/N). The normalization 
for the finite- A" corrected calculation deviates from 1 within ±1 per cent. 

Figure 2 presents detailed comparison between Monte-Carlo results and the above 
analytical approximations for three different values of E av , {—0.95, 0, 0.5}. The numbers 
of Monte Carlo points accepted by the algorithm from Appendix B for the three values of 
E av were, respectively, 2000, 5000 and 5000. The first row of frames in Fig. 2 represents 
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Figure 2. Comparison between the results of direct Monte-Carlo sampling of the 
QMC ensemble for N — 12 quantum spectrum given in the text and in Fig. 1 (blue 
circles with statistical error bars, when visible) and three analytical approximations 
based on functions A[-E a v] presented in Fig. 1: magenta triangles and magenta long- 
dash lines — small-pfe approximation based on Eqs. (|15|16|) : green diamonds and green 
small-dash line — mean- A approximation based on Eq.(|27j): red squares and red solid 
lines — finite- N corrected calculation based on Eq. ([26| and the ansatz of Section [5l 
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the comparison for the values of the average occupation numbers (pk) corresponding to 
twelve energies in the spectrum. The remaining three rows represent the comparison 
for the three (out of twelve) probability distributions Pk(pk) for each of the above three 
average energies. 

In the case of P av = —0.95, the Monte-Carlo result qualitatively corresponds to 
the prediction of Ref.pQ that, for macroscopic systems (meaning Gaussian density of 
states in the large- N limit) at realistic positive Hilbert space temperatures 1/(NX) (i.e. 
with P av < E av0 ), the QMC ensemble leads to condensation into the lowest energy 
state with P\(pi) having the character of ^-function peaked at (pi), while all other 
states have small average occupations (pk) and exhibit exponential shapes of Pk{pk)- 
What one sees in Figs. 2(a,b,c,d) is a finite- N realization of the above condensation. 
The Monte-Carlo results indicate that all Pkipk) except for Pi(pi) corresponding to 
the lowest spectral level E\ = —1.27 have nearly exponential shape well describable 
by all three of the approximations considered [Fig. 2(c,d)]. The differences, however, 
become pronounced once Pi(pi) is considered. Figure 2(b) indicates, that the small-p*,. 
approximation fails completely to describe Pi(pi), the mean A approximation reproduces 
the Monte-Carlo results qualitatively, but with large quantitative discrepancy, while the 
finite- N corrected calculation is in a good quantitative agreement with the Monte-Carlo 
results. 

Despite failing to predict Pi(pi), the small-p^ approximation predicts accurately 
the value of (pi) in Fig 2(a). As discussed in Ref. [TJ, it is the consequence of the fact 
that this approximation accurately describes all Pk(pk) and (pk) for k > 1, but then, 
since it automatically satisfies the normalization condition ( ITT]) , it guarantees that the 
value of (pi) is also accurate. 

The example of P av = 0.5 corresponds to the negative Hilbert space temperature 
l/(iVA). In this case, the QMC ensemble tends to condense to the highest energy 
state. The facts, that, in this case, the distribution of ^12(^12) [Fig. 2(1)] is rather 
broad, and that the form of Pn(pn) [Fig. 2(k)] for the second highest level also 
noticeably deviates from exponential, are the finite-iV effects. Here again, the small- 
Pk approximation fails to describe Pn(pn) and Pi2(pi2) , the mean-A calculation gives 
correct qualitative trends with quantitative errors, and the finite- iV-corrected calculation 
exhibits a good quantitative agreement with the Monte-Carlo results. For k < 10 all 
three approximations give nearly exponential decays for Pk{Pk)i i n good agreement with 
the Monte-Carlo results [see, e.g., Fig. 2(j)]. Since the small-p^ approximation fails for 
three highest levels, it is no longer guaranteed to produce the accurate values of (pio), 
(pn) and (P12) and, in fact, the noticeable discrepancy can be seen in Fig. 2(i). 

The example of P av = corresponds to the infinite Hilbert space temperature. In 
this case, the small-p^ approximation predicts that all Pk{Pk) have the same exponential 
shape, and that, independently of (pk) = l/N. This would indeed be the case in 
the limit iV — > 00. However, here the symmetric position of P av = with respect to the 
quantum spectrum leads to a more systematic character of finite- N corrections, which, 
in turn, leads to a larger discrepancy with the Monte-Carlo results for all k. The origin 
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of this systematic discrepancy is in the corrections to the linearized version of integral 
in Eq. ( |2"7|) . The linearization of that integral with respect of pk leads to the small-p^ 
approximation. The higher order corrections are, in general, partly canceled, because 
X[E] in the integral changes sign, when the integration range includes E = somewhere 
in the middle. However, when E av = 0, this sign-change point coincides with the lower 
integration limit. As a result, X[E] does not change sign within the integration range for 
any k, which leads to the above larger discrepancy. Two other approximations capture 
the above systematic trend, but here again, the finite- N corrected calculation exhibits 
a better quantitative agreement with the Monte-Carlo results. 

It should be emphasized that the finite-iV corrected calculation is based on an 
ansatz of Section [5], which, although meaningful and apparently very effective, does not 
amount to an exact calculation. Therefore, it is expected to show small deviations from 
the Monte-Carlo results, which can, indeed, be seen in Figs.2(k,l). 

The statistics of eigenstate participation for the entire isolated system reported in 
this section is not a direct test of the presence or the absence of BG equilibrium for a 
small subsystem. However, in this respect, the results of significance are those shown 
in Figs. 2(a,e,i). They represent, in a sense, the participation functions for conventional 
micro-canonical ensembles with different positions of micro-canonical energy windows. 
These participation functions have algebraic dependence on energy mostly captured by 
the small-pfc result (fl6l) . This slow non-exponential dependence precludes one from 
treating the participation of eigenstates as a conventional canonical ensemble. As a 
result, the broad participation of eigenstates in the QMC ensemble implies a broad 
mixture of conventional thermal states, which, for a small subsystem, does not result in 
a single thermal BG distribution — the conclusion confirmed by direct calculation in 
Ref.pp. In the next section, we present numerical evidence corroborating this conclusion. 

7. Implications of the QMC ensemble for a subsystem within a small spin 
system 

Reference [1] derived analytical results for the density matrix of a small subsystem 
within a larger isolated quantum system described by the QMC ensemble. When the 
whole system is macroscopic (i.e. the number of levels N is exponentially large), the 
QMC-based density matrix of the small subsystem is a weighed sum of two terms: 
the zero-temperature density matrix and the infinite temperature density matrix. In 
the finite-iV case, one has to solve a system of equations obtained in pQ ■ - not to be 
attempted in the present work. Instead, in this section, we demonstrate numerically that 
the following qualitative property of the macroscopic limit also appears in the finite- N 
case. Namely, for a subsystem within an isolated quantum system, the lowest energy 
state and the states on the high-energy end of the spectrum have larger occupations when 
the isolated system is sampled according to the QMC ensemble than when it is sampled 
according to the conventional canonical ensemble. Correspondingly, the intermediate 
energy states have smaller occupations for the QMC ensemble. 
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Figure 3. Energy spectrum of the 12 level system (two spins 1/2 plus spin 1) 
considered in Section [7J 

Below we consider a 12-level system consisting of two spins 1/2 and spin 1, and 
examine numerically the implications of the QMC ensemble for the subsystem of two 
spins 1/2. 

As discussed in the previous section, the 12-level limitation is due to the 
computational difficulty of sampling high-dimensional Hilbert spaces. In principle, in 
order to exhibit deviations from the BG statistics, the subsystem should have at least 
three levels, and therefore one could think that it is more natural to consider the spin 
1 as a subsystem. However, with only 3 levels in the subsystem the differences between 
the QMC-based results and the BG statistics would be less instructive. 

We assume the the environment represented by spin 1 does not interact with the 
subsystem of two coupled spins 1/2. Hence the Hamiltonian of the entire system is: 

U = U I + U S , (39) 

where 

Uj = nH{I lz + I 2z ) + JJ lxhx + Jyhyhy + JJuhz (40) 

and 

Us = isHS z (41) 

are, respectively, the Hamiltonians for the two spins 1/2 and for spin 1. Here I ma are the 
operators of for the a-components (a = x, y, z) of the two spins 1/2 (m = 1, 2), S z is the 
operator of the z-component of spin 1, H — 3 is the external magnetic field, 7/ = 0.7 and 
7s = 3 are the gyromagnetic ratios for spins 1/2 and for spin 1, respectively; J x = —2, 
J y = — 1, J z = 0.5 are anisotropic coupling constants. Since spin 1 plays the role of the 
environment, the parameters of the Hamiltonian were chosen such that the spread of 
the spin-1 levels is significantly broader than that of the two-spin-1/2 subsystem. The 
energy spectrum of the Hamiltonian fl39|) with the above choice of parameters is shown 
in Fig. [3j It has three clearly identifiable 4-level groups. 

We compare the subsystem density matrices for the QMC ensemble and for the 
canonical ensemble at the same value of the average energy of the entire isolated system. 
In order to do this, we first fixed the temperature of the canonical ensemble T=l, 
then calculated the corresponding average energy E m = —10.55, and then generated 
2000 sample wave functions from the QMC ensemble for the above value of -Eav- 
Afterwards, we calculated the average subsystem density matrices corresponding to the 
both ensembles. The results for the diagonal elements of the density matrix p aa in the 
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Figure 4. Diagonal elements of ensemble-averaged density matrix < p aa > for the 
subsystem of two spins 1/2 considered in Section [7J in the basis of eigenstates of the 
Hamiltonian Hi. The subsystem eigenstates are labeled by index a with corresponding 
energies denoted as E a . Magenta squares represent the QMC ensemble, blue circles 
represent the canonical ensemble. Inset: same results on the linear scale. 



eigenbasis of the Hamiltonian Hi are presented in Fig. HI As expected from the large- iV 
limit p], the values of p aa for the QMC ensemble are larger than those for the canonical 
ensemble for the lowest and the highest subsystem energy levels, and correspondingly 
smaller for the middle two levels. 

If one can have only limited experimental access to the properties of the subsystem, 
the QMC ensemble can be identified by looking at the variance of energies for the 
subsystem states: 

(AE P ) 2 = J2p^(E a - E p )\ (42) 

a 

where index a labels subsystem's eigenstates, E a denotes the corresponding subsystem 
eigenenergies, and E p = J2 a p aa E a is the average energy of the subsystem. In Fig. [5j we 
plot subsystem energy variances for the QMC and the canonical ensembles as a function 
of E p . For E a y = E\ (see Fig. |3]) and E m = 0, the predictions of the two ensembles 
are identical, and hence the variances are the same. (In the first case, only the ground 
state is occupied and E p = min[-E a ], while in the second case all subsystem's states are 
equally occupied and E p — 0.) In-between the energy variance for the QMC ensemble 
is always larger. This is true not only for the subsystem, but also for the entire system. 

An interesting fact is that for the parameters considered, the subsystem average 
energies for the QMC and the canonical ensembles are different even though the energy 
of the entire system is the same for the both ensembles. The subsystem energies for the 
two ensembles are, respectively, E® MC = —1.83 and E c p ^ n = —1.55, i.e., in the QMC 
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Figure 5. Energy variances for the subsystem of two spins 1/2 given by Eq. (|42|) as 
a function of the average subsystem's energy. Magenta squares represent the QMC 
ensemble, blue circles represent the canonical ensemble. 



case, the subsystem is in a sense "colder". We have investigated this issue further, and 
in Fig. [6] have plotted the result for E p for the both ensembles as a function of E av of 
the total system. Our findings indicate that there exists a critical value E™ ~ —4 such 
that E^ MC < E™ n for E w < E™, and E® MC > E c p an for < £ av < 0. However for 
< E av < 0, the results for the two ensembles are very close (the high-temperature 
limit) and hence the difference between E® MC and E^ an is insignificant and difficult to 
resolve statistically. 

The difference between E® MC and E c ^ n may, perhaps, be used as another 
experimentally identifiable signature of the QMC ensemble. It should be noted, however, 
that the above difference can only appear when the subsystem particles are different from 
the environment particles. The subsystem energies for the two ensembles are guaranteed 
to be the same, when, for example, the whole system is describable by a translationally 
invariant Hamiltonian of spins 1/2, and the interaction between the subsystem of a few 
spins 1/2 and the environment is negligible. 

8. Prospects 

Here we review relevant open questions. 
8.1. Macroscopic systems 

It appears now very likely that the micro-canonical narrow energy window condition 
is crucial for the justification of BG equilibrium in macroscopic systems. What, 
however, can justify this condition itself? The only argument known to us is associated 
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Figure 6. Average energy for the subsystem of two spins 1/2 as a function of the 
average energy for the whole system. Magenta squares represent the QMC ensemble, 
blue circles represent the canonical ensemble. Inset: Ratio of the subsystem energy for 
the QMC over the subsystem energy for the canonical ensemble. 



with quantum collapse - - also referred to as quantum measurement or quantum 
projection[l] |8]. Eigenstates with significantly different energy values should have 
macroscopically distinguishable characteristics. Hence the mixture of macroscopically 
distinguishable states, coherent or not, would not be tolerated in nature, and instead 
the mixture would collapse to one macroscopic state with usual probabilities predicted 
by quantum mechanics. 

The disturbing part of the above argument is that the quantum collapse is an 
essentially non-linear process, and hence cannot appear in the quantum-mechanical 
theory in the course of linear evolution of a wave function or a density matrix. Therefore, 
it implies that BG equilibrium cannot be justified in the framework of linear quantum 
mechanics — it requires a non-linear quantum collapse into the narrow micro-canonical 
energy window. 

One can still ask a different question: If the quantum collapse occurred only 
once in the evolution of a macroscopic system, is it realistic to expect that external 
perturbations on otherwise isolated macroscopic quantum system would force the spread 
of participating eigenstates over an energy window violating condition (J2])? Is it further 
realistic to expect that the QMC ensemble can be realized in an isolated macroscopic 
quantum system that has once experienced quantum collapse? For a truly macroscopic 
system consisting of 10 23 weakly interacting parts, the answer to both questions appears 
to be "extremely unlikely". The reasoning here is based on the central limit theorem: 
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If such a system is subjected even to a drastic external perturbation, e.g. explosion, 
it is still likely that different parts of such a system, while being strongly perturbed, 
will continue to interact weakly, end hence the resulting total energy distribution will 
have Gaussian form, which would be sufficiently narrow to satisfy condition (j2]).(See 
also Refs.pHHSl E].) 

Yet a single perturbation makes the energy window of participating eigenstates 
somewhat broader. The next level of this argument would be to ask, how many 
sequential perturbations are necessary to force an isolated macroscopic system to violate 
condition <^j. This question requires further investigation. Realistic macroscopic 
systems are continuously subjected to the fluctuations of external potentials of 
electromagnetic or gravitational origin. Can these fluctuations drive an isolated system 
out of the micro-canonical energy window? If they can, but we do not observe any 
manifestation of it, does it mean that quantum collapse to the narrow energy window 
occurs continuously? If this happens, does it imply observable energy fluctuations not 
explainable by linear quantum mechanics? 

8.2. Small quantum systems with many levels 

The second line of investigation is related to dealing with quantum systems that do not 
yet have large number of degrees of freedom, but already have a very large number of 
quantum levels, e.g. 10 x 10 x 10 clusters of interacting spins 1/2 or a few cold atoms 
in an optical trap. 

In this case, it appears realistic to isolate these systems against energy dissipation 
and, possibly, against decoherence. The central-limit-theorem-based considerations are 
not applicable to these systems. Therefore, it is clear that such systems can be driven out 
of condition described by conventional micro-canonical or canonical ensembles. Will the 
QMC ensemble generically emerge in this case, possibly, after several perturbations [16]? 
If such an ensemble emerges, what are the observable signatures of it, e.g. in terms 
of measurable single particle properties? What would be the effects of decoherence 
and subsequent quantum collapse in this case? It is, in particular, possible that even 
if the broad QMC energy window collapses towards the canonical or micro-canonical 
shape, then the QMC-based calculations would predict the fluctuations of the resulting 
temperature. These questions require numerical investigations beyond the scope of this 
work. It is, in fact, likely that non-micro-canonical ensembles, and possibly QMC, 
routinely occur in the numerical studies of finite quantum systems. 

9. Conclusions 

In conclusion, we have demonstrated that the large-iV analytical description of the 
QMC ensemble developed in Ref . [1] allows one to make good predictions for the Monte- 
Carlo sampling of the QMC ensemble for finite- iV systems. The description of Ref. pQ is 
amenable to finite- N corrections, and once these corrections are introduced, the theory 
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produces accurate quantitative agreement with the Monte-Carlo results for systems with 
N ~ 10. 

We have also studied numerically the implications of the QMC ensemble for a 
subsystem of a small spin system. Our results indicate that, already in this rather 
artificial case, the behavior of subsystem's density matrix qualitatively follows the 
analytical result for the large-A case. We have further suggested that the subsystem 
energy variance, as well as the average value of subsystem's energy itself, can be used 
to discriminate experimentally between the QMC and the BG statistics. 

In the context of foundations of quantum statistical physics, the apparent 
impossibility of obtaining the Boltzmann-Gibbs equilibrium from the QMC ensemble 
indicates that one should be looking closer at the basis and the applicability limit of the 
narrow energy window condition of the conventional micro-canonical ensemble. In this 
respect, quantum collapse may play a fundamental role. 

Note added: As we were finishing editing this manuscript for resubmission, we 
discovered a series of papers by Fresch and Moro[TFl HHJ [19] very closely related to the 
scope of the present work and Ref.pQ. In particular, in Ref.[T7], the above authors 
have pursued similar Monte-Carlo sampling of the QMC ensemble in finite dimensional 
Hilbert spaces. Their results appear to be consistent with and, in some aspects, 
complementary to ours. In particular, they have used Metropoilis-Hastings algorithm, as 
opposed to the direct sampling routine of our paper, and thereby were able to explore 
significantly larger Hilbert spaces. On the other hand, the present paper advances 
further the subject of finite-A corrections to the large-A QMC results. We have also 
learned from Ref.[T7] that, already in 1990, Wootters have considered the QMC ensemble 
and arrived to what we call the small-p^ approximation [20] . 

Appendix A. Algorithms for finding A[i? av ] 

Here we describe the algorithm for finding A[-E av ] from the mean-A integral equation ( |27j) 
- "Algorithm A" , and the modification of Algorithm A including finite- A corrections 
based on Eq. (l26l) incorporating ansatz from Section |5] — "Algorithm B". 

Function A[i? av ] is defined in the interval [E min , E max \. We use a non- uniform 
discretization of this interval with 228 grid points - - less dense in the middle and 
more dense near E m[n and E max , were A[-E av ] tends to diverge. This near-divergence is, 
however, manageable, because, it is compensated by the overall exponential decrease of 
Vkipk) and affects only the range of pk, where Vkipk) is exponentially small. 

Algorithm A: Outline 

The algorithm first finds A[-£7 av ] in the small-p^ approximation on the basis 
of Eqs.(fT6l) and ( fi~8l) and then iteratively improves it by (i) calculating Vk(pk) by 
substituting A[i? av ] into Eq. fT27|) . (ii) finding the corresponding (pk), (hi) checking 
averaged energy condition ( JTBl ; (iv) improving the value A[i? av ] by matching condition 
f fT8|) and then returning to step (i). As explained below, step (iv) is done "locally", i.e., 
for every E av , on the grid it improves A[-E av ] independently of the values of A[£7 av ] at 
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other grid points, but then once all grid values of A[i? av ] are improved, the algorithm 
returns to step (i) with the entire improved set of A[£" av ]. 
Algorithm A: Details and explanations 

At the small-pfc approximation step, we substitute {p k ) from Eq. (fl6l) into Eq. (fl8]) . 
and then find numerically the value of A such that i? av — j- < E min or E av — j > E max . 

At the first iteration step, we start by substituting A[-B av ] obtained in the small- 
Pk approximation into the integral formula (j2"7|) . thereby obtaining the generally non- 
exponential function V k (p k ), then use Eq. ffTTj) to obtain (p k ) and then check condition 
(Fl8|) . which shows deviation from zero. In response to this deviation, we modify the 
function A[i£ av ]. 

The value of A[i? av ] at one grid point affects, according Eq. fl27|) . the values of V k {p k ) 
and hence condition ( fl8l) at other grid points. Therefore, the proper iteration requires 
finding the whole matrix of linear responses of function 

N 

F(E av ) = J2(E k -E^)(p k ) (A.l) 

fc=i 

at each grid point to the change of A at each other grid point. If the _E av -grid has Nq 
points, then one has to solve a system of Nq linear equations to find the modification 
AA[i? av ] for each grid point. From our experience such an approach often produces 
spurious solutions and, therefore, is difficult to manage. Instead, we adopt a different 
approach based on our understanding, why the small-p^ approximation works so well. 

In the small-pfc limit, the algorithm for finding A[_E a v] is "local", i.e. it does not 
involve explicit influence of A[i2 av ] at a given grid point on F(E' SV ) at other grid points. 
In this case, the value of A is tried, and if it does not give F(E av ) = 0, then it is modified 
to A' = A + AA, which, according to Eq.f JT5|) . implies that V k {p k ) is also modified to 
become 

VLiPk) = Vk (0)e- Np ^ x+Ax ^- E ^. (A.2) 

This equation can now be rewritten as 

VM = V k ( Pk )e- NpkAx ^- E ^\ (A.3) 

In the above form, the iteration step is generalizable to the case, when V k (p k ) itself 
deviates from the exponential form. We find the value of AA such that F(E av ) calculated 
with (p k ) corresponding to Vl(p k ) substituted into Eq. tflTT) equals zero. Such an ansatz 
for V k (p k ) is not fully consistent with Eq. fT27|) . but it is a good approximation of it, and 
more importantly leads to a good improvement in the value of A. The version of V k (p k ) 
consistent with Eq.( )27|) is then obtained at the beginning of the next iteration step 
by substituting the entire improved function A[i? av ] into Eq.(|271). which, in turn, leads 
to the deviation of F(E ay ) from zero, which is again canceled at each E av by replacing 
Vkipk) by V k {p k ) with the new value A A, etc. In our experience, the form of A[i? av ] stops 
noticeably changing after three such iteration cycles. The calculated results presented 
in Section |6] are based on six iterations. Ultimately, it is the observed quick convergence 
that justifies the validity of the present algorithm. 
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The rationale for updating the values of A[J5 av ] on the basis of Eq. (1A.3|) is that 
(i) change of A at one value of E w affects mostly F(E' av ) for E' av reasonably close to 
E^; and (ii) for E' av close to E av , the value of AA[£7 av ] should not be very different 
from AAfE^v] - hence AA[_E av ] well approximates AA[-E^ V ] once the corrections are 
substituted in Eq. (l271) . The correction AX[E av ] obtained with the help of Eq.( 1A.3jl 
accurately anticipates the linear term in the expansion of integral in Eq. (l27[) in powers 
of pk and the resulting exponential change of Vk{Pk) for sufficiently small values of pk- 
For the larger values of pk, the exponential factor in Eq.l 1A.3l) introduces a larger relative 
error, but, in almost all cases, this is a larger error to the small and hence less important 
tail of Vk(pk)- The tail is then corrected at the beginning of the next iteration step once 
the updated function A[-Ea V ] is substituted into Eq. 11271) . 

Algorithm B 

Algorithm A can now be modified to incorporate the finite- A" corrections associated 
with the deviations of individual functions Afc[i£ av ] from the "mean" function A [E'av]. 

Algorithm B continues searching for mean function A [E'av], but at each step of 
Algorithm A that used Eq. (l2"T|) . it uses Eq.f l26|) with values of X^E^] obtained from 
A [-E'av] with the help of ansatz described in Sectional The algorithm continues using the 
"local" step involving Eq. (1A.3j) . and, here again, the fast convergence of the iteration 
procedure is the evidence of the adequacy of that step. 

Appendix B. Algorithm for direct Monte-Carlo sampling in 
energy-constrained Hilbert space 

Without energy constraint ([6]), there exist a straightforward and efficient algorithm 
for the sampling of Hilbert space. Namely, one can first choose all values of Re(Cj) 
and Im(Cj) according to the same Gaussian probability distribution. This gives 
unnormalized wave function. After it is normalized, the resulting Monte-Carlo points 
have the required isotropic probability density on the normalization manifold (provided 
one can control the effects of the cutoff on the tail of Gaussian distribution necessarily 
appearing in the numerical implementation of the sampling). 

With energy constraint ([6]), we are not aware of any algorithm, which would be 
as efficient as the above one and yet rigorously implement uniform sampling on the 
energy-normalization intersection manifold. 

We used the same algorithm as Ref.[T]. Namely, we do the sampling in the {pi}- 
Euclidean space, where we define an orthogonal basis in the (N—2) dimensional subspace 
formed by the intersection of energy and normalization hyperplanes ([6]) and (J7|). The 
orientation of these hyperplanes are determined by the unit vectors, which are orthogonal 
to them: w E = * (Ex, E 2 , E N ) and v norm = hL(1, 1, 1). The basis 

within the intersection manifold is then generated by Gramm-Schmidt orthogonalization 
procedure. Namely, we make the first basis vector v x orthogonal to both v E and v norm ; 
then v 2 orthogonal to v E , v norm and Vj.; etc. Each of (N — 2) vectors Vj then gives rise 
to a new coordinate axis, which we denote as q±. One has, of course, much freedom in 
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choosing this basis. 

In the above (N — 2)-dimensional subspace, the QMC polyhedron is further 
constrained by the positivity conditions pi > 0. As a result, it has KL vertices [T], 
where K, and L are the numbers of energy levels in spectrum {Ei} above and below 
E av , respectively In the original basis of space {pi}, each vertex has only two non- 
zero coordinates (0, 0, p m , 0, 0, p n , 0, 0) corresponding to one possible pair of 
energy levels such that E m < E av and E n > E av . The values of these coordinates are 



We choose the origin of the (^-coordinate system on one of the above vertices. The 
Pi coordinates of each vertex are then transformed into qi coordinates. Then, in the q- 
coordinate system, we choose an (N — 2)-dimensional hyperrectangle (orthotope) with 
edges parallel to the qi axes and with the extension along each axis limited by qWin and 
9imax such that gj-coordinate of each polyhedron vertex falls between gj min and gj max , i-e. 
the polyhedron of interest is completely inside the above hyperrectangle. 

The Monte-Carlo sampling is then straightforwardly implemented in the 
hyperrectangle by choosing each coordinate qi randomly in the interval [q^mm, ^J. The 
Monte-Carlo point is accepted, when, after transformation to the original pj-coordinate 
system, all of its coordinates satisfy the positivity condition p { > 0. 

The acceptance rate of this algorithm decreases exponentially with increasing N, 
because high-dimensional polyhedron typically occupies exponentially small part of the 
volume of a hyperrectangle that covers it. 

The above ratio of volumes has a further strong dependence on the relative 
orientation of the hyperrectangle and the QMC polyhedron. It was our guess confirmed 
empirically that this ratio is strongly increased, once one of the faces of the polyhedron 
is parallel to one of the faces of the hyperrectangle. The faces of the polyhedron 
are determined in the p^-space by the intersection of three hyperplanes: Pi = 0, 
energy hyperplane (EJ) and normalization hyperplane (CO). If one of vectors v« is 
orthogonal to one of the polyhedron faces, then that face will be parallel to the 
corresponding face of the hyperrectangle. For example, if we want vi to be orthogonal 
to the face determined by hyperplane p\ — 0, then we choose v x to be parallel to 
ui - (ui • v E )v E - (ui • v norm )v norm , where u x = (1, 0, 0, 0) in the original p^-space. 
For N = 12, the resulting increase in algorithm's acceptance rate was, approximately, 
by a factor of 50. 

We have not yet systematically optimized the choice of the "parallel" polyhedron 
face as a function of the corresponding energy E^ from the viewpoint of increasing 
the polyhedron-to-hyperrectangle volume ratio, but, so far, our experience indicates 
that, for E av < E av0 , the face corresponding to the lowest energy, i.e. originating from 
hyperplane p\ = 0, is one of the best. The algorithm can, perhaps, be further improved 
by optimizing the subsequent choice of V2, V3, etc. from the viewpoint of increasing the 
above volume ratio. 



Pm = 



E n -E m alld ^ 
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